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UNIVERSITY OF ILLINOIS 
DIGITAL COMPUTER 



LIBRARY ROUTINE L 8 - 302 



Solution of a system of linear equations by an iterative 

method (SADOI Only) 

Complete program 

338 (20-358) in Williams Memory 

Depends on the condition of the equations. For more 

information on this point see the description of the 

mathematical method. 

Depends to a large extent on the initial approximation. 

If it's close to the exact solution, convergence will be 

fairly rapid; if not, then convergence is liable to be 

very slow. If T is the time required to solve a system 

of N equations in N unknown, then T - .013TT minutes. For 

example, when N = 35, T = 16. This estimate of running time 

is to be regarded as an approximation only! 

S6" - OOF OOnF where n = number of equations 

S9 - OOF OOkF where 1 < k < 12 is the number of 

decimal places desired in the output* 
A must be symmetric and positive definite. In addition 



N-l 

Z |a 
j=0 



ij 



< 1 for i = 0, 



N-l. 



This implies that the largest eigenvalue of A < 1. On the 
data tape, to be described below, the user will have to 
specify a and b where a is the lower bound for the magnitude 
of the eigenvalues of A and b is an upper bound. cannot 
be taken as a lower bound for the eigenvalues of A. Finally 

N < 100. 

If there has been a drum failure there will be an FF stop 
at R.H. (0F0) r otherwise, when the problem is done, there 
will be an OF stop at (0S9 )-[_£. 



- 2 - 



METHOD OF USE: 



Suppose the set of equations to "be solved is A x = y where 
A is an M x M positive definite symmetric matrix and y a 
vector with M components. Let a denote the magnitude of 
the smallest eigenvalue or a lover bound for the magnitude 
of the smallest eigenvalue and let b denote the magnitude of 
the largest eigenvalue or an upper bound for the magnitude of 
the largest eigenvalue of A. Library routines M 26 or M 20 
may be used to determine these bounds. 
It is essential that 



< a c b < 1. 

The data tape is then prepared as follows: 
a 
b 
s this is a scaling factor, to be discussed below 



11 
*12 



l uJ 



.rst row of matrix, 



each row of matrix must 



21 

So 



2Mj 



be followed by 

the sexadecimal character N, 



a "\ 

Ml 






Last row of matrix, 



Followed by: 



s o y i 



SCALING: 



wnere y 



(y l '•• ^ 



8 y MJ 

N 



S 0^01 



S 0?0M 
N 



Followed finally by: 



where h . 







- (1 



01- 



;0M 



) is the 



initial approximation 



All the above numbers are signed decimal fractions which 

are read in via N 12 . 

-1 -2 
s« will usually be 10 or 10 etc. but other scale 




factors such as 2 



-1 



-2 



etc. may also be used. The 



scale factor should be chosen so that the matrix s_ A 

satisfies the restrictions noted above and the coordinate 

elements of s y and s C Q satisfy the following inequalities 

(in order to prevent overflow): 
M 

< 1/2 



1) 



s o Z 1 



'0j 



2) 



M 



< 1/2 



The routine automatically performs any additional scaling. 
The first number punched on the output tape is s, _, the 
final scaling factor, followed by s, x, the scaled solution 
vector. 
BRIEF DESCRIPTION OF MATHEMATICAL METHOD AND CONVERGENCE CRITERION: 

Denote the system of equation to be solved by 

(1) A? = ? 



- k - 



It Is assumed throughout the following discussion that 
A is symmetric, positive definite, and has all eigenvalues 
less than 1. Now if one wants to solve this system of 
equations by an iterative method then, as has been pointed 
out "by von Neumann [2], a large number of such schemes are 
described by the equations: 



< 2) f w = S ?* + v 



where | denotes the kth iterate and G, and.H, are matrices 

of the same order as A. The G and K are not independent 

of one another; they must satisfy the following two conditions: 

(3) G k + H k A =1 (I - identity matrix) 
00 det iHj £0 



The reason for this is the following: suppose 2 is 
the solution i.e. \ = x then obviously "5 = y 
which means x = (G, + KA) x, all x, hence (3) follows. 
Similar considerations lead to condition 4. For further 
details consult Golub [l] pp. 7-8. Since x = (G + H,A) x 
then upon subtracting (2): 

(6) *" fk+1 ssG k G? " Tk } (since A? =y*). 

Therefore 

(7) (?-1 k ) -*" 1 G,(?- f„). 

> k j-0 J 5 ° 

In many cases G, = G(all k) so (7) becomes 

(8) ?- ^ k = aV-j ). 



_ ■=> - 



It is now easily seen that if the eigenvalues of G are 



-> 



all less than 1 in magnitude x - j? -» independently 
of ? n . However, if ^ is a poor approximation or if 
the largest eigenvalue of G is close to 1 the convergence 
may be too slow. In [2] von Neumann considers replacing 
the original sequence ^ , -^ , by a sequence of averages 
~n > "?p ••• vn ere 



k 



&)■ fk% Z 



k £=o k 



a, * 



J )L 



and 

k 
(10) ,E a vl? = 1. 
=0 



i=e~ kiL 



Hence 

k. =0 ' -*-' / =0 



L a ^ 



=o •" '° 



That is, after k iterations of this process the initial 

error is multiplied by a matrix polynomial so that 

k n 
(12) x - "I = P k (G) (x - f ) where P k (G) = £ \p^' 

jC- =0 

If. x - ~$ _■ is expanded in terms of the eigenvectors of 
c 

G(if G is symmetric, this can always be done) then the 
following inequality results: 



(13) ||x- ? ||<max IMMI ' I I* " f 

k ~ 1< i < N K X 



6 - 



where k. are the eignevalues of G. Now for arbitrary ~$ 



| jx - h' | | is to be as small as possible after k steps, 
hence max |P, (X. )| must be minimized. The answer is well 
known and leads to the Chebyschev iteration scheme described 
by Golub [l] pp. 9-15 • The general expression isi 

(Ik) ^ +1 . 2b k+1 (G .? k + tf - f^) + ^ ( k . 1,2, ...)■ 

b k+l = — ■ 2 — » ^t - 1 and - *• is an u PPer bound for 
2-X, b. 
■ k 

the largest eigenvalue . ..._'-■■ 

In the library routine described here H. = nCl, all k, and 

hence G =1 - ^ A, where o4 is chosen so as to minimize 

the maximum eigenvalue of G. If the smallest and the 

largest eigenvalues of A are denoted by a and b respectively, 

then *S = — £. This choice of o^ yields the iterative method 

used in this routine. 

CONVERGENCE CRITERION: 

As is pointed out by Golub, | \to - h \\ < e does not 

imply | |x - *7fc + ill < e * However, if z! .is defined as 

\ = G fk + H ^" fk = m (* " ?k^ then if I l ? k' I is 
"small" so is | |x - K? k | | • More precisely 

-lii 



- - Vj k so ||x- ^J| < 



-T A z k = x - ^ k so | |x - TJ k | | < ^-u. | | Zk 



Therefore if 



--> ■ - - JaLL 



k " MA" 1 



- 7 - 
then 

II?- ^ k H< e - 

Now 
EA 



6?. ^ k )«o<A<?- f k ) =<*(?- A ? k ) = *? k = V 

Because of round off errors e cannot be made arbitrarily 

small. In fact it might be the case that min| |z, | | > cf > 

k 

for all k. Hence, if € is chosen too small, the process 

may not converge. A lower bound for € has been determined 

by Golub (see [l] P. 86) which can be expressed as 

M 

€ < 



d-cy 

where 

e = - " 



M = TO • 2 Jy 



and 

b-a 



X = fc+ft / 

M is an estimate of the rounding error which occurs during 
iteration. It should be noted that if C is close to 1 
(which is the case for ill conditioned matrices) e becomes 
large. Therefore, much less accuracy is to be expected 
in this case. 



- b 



*i£F 



;ss, 



[1] Digital Computer Laboratory, University of Illinois 
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LOCATION 




ORDER 




NOTES PAGE 1 L 8 






00 1QK 












26 1000N 













L5 6F 
00 IF 




36l + 2N 




1 




1^ 10L 
k2 3F 




(S3) 




2 




50 6f 
75 6f 








3 




L5 6f 
S4 P 




N 2 + N + 2560 




k 




lA 11L 
40 ^F 




(s h) 




5 




L5 6f 
00 IF 




3N + 1 + 360 




6 




1,4 6f 
1^ 10L 








7 




^2 5F 
L5 12L 




(S3) 




8 




k2 JF 
L5 13L 








9 




42 8F 
26 999F 








10 




00 F 
00 56lF 








11 




00 F 

00 2560F 








12 




00 F 
00 202F 








15 




00 F 

00 251F 

26 ION 
00 20K 









r 



LOCATION 



3 
4 

5 
6 

7 
8 

9 
10 
11 
12 

13 
14 

15 
16 



(o) 



ORDER 
00K(L9) 

40 (506) 

L5 OF 
10 IF 
40 OF 
L5 IF 
10 IF 
40 IF 
L4 OF 

4o (507) 

L5 IF 
LO OF 
50 (508) 
66 (507) 
S5 F 

4o (509) 
75 (509) 

10 IF 

40 (510) 
L5 (501) 
L4 (505) 
40 (503) 
L5 (500) 
L4 (505) 
40 (502) 
L4 (511) 
40 (504) 
50 (509) 
71 (509) 
L4 (511) 
40 IF 
S5 F 
40 OF 
00 IF 
50 (0) 



NOTES 
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Store scale parameter 



i/2 



b/2 

a/2 + b/2 

Clear Q 

\ = b-a/b+a 

\ 2 /2 



| LOCATION ' 




ORDER 




NOTES PAGE 3 


L8| 


17 




22 (Rl) 
10 IF 




n - y? 




18 




40 OF 
LT OF 








19 




40 OF 
50 (508) 








20 




L5 (509) 

10 IF 








21 




66 of 
35 F 




<° = \/i + iZ>F 




22 




L4 (511) 

4o OF 








25 




50 OF 
7J OF 




d-ef 




24 




40 OF 
L5 (512) 








25 




50 (508) 

66 of 








26 




ST F 
40 (513) 




Compute end constant 




27 


(25) 


49 OF 
LI OF 


from 108 






28 




40 (514) 

41 (523) 








29 


(2h) 


22 (1) 
L5 (514) 


from 90, 99 






50 




00 IF 
36 (2) 




Is -b = -l/2? 




31 




89 IF 
40 (514) 








32 


(2) 


50 (514) 
79 (510) 


from 30 






33 




L4 (511) 












40 OF 








34 




50 (508) 












LJ (508) 









LOCATION 
35 

36 

37 

38 

39 

1*0 

41 

42 

43 
44 

45 

46 

47 
48 
49 
50 
51 
52 



(1) 



(17) 



(3) 
(7) 



(18) 



ORDER 
66 OF 
S5 F 

40 (514) 

41 (517) 
41 (515) 

41 (516) 
L5 (500) 
40 (3) 
40 (19) 
L5 (501) 
40 (4) 
L5 (518) 
40 (5) 

42 (6) 
L5 (520) 
42 (7) 
42 (8) 
42 (9) 
46 (6) 
42 (11) 
85 11F 
40 S4 
32 (7) 
40 S5 
F5 (3) 
40 (3) 
LO (502) 
36 (18) 
F5 (7) 

te.(7) 
LO (521) 

36 (3) 

L5 (6) 
40 IF 

L5 (517) 
40 OF 



NOTES 



PAGE 4 



L8 



from 29 



Compute - b 



k+1 



from 82 



"by 38,47; from 5° 



by 42,49 



Read in 128 components of A h , ^ 



from 48,75 



- 

LOCATION 




ORDER 
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55 


(10) 


L5 (506) 
50 (10) 






54 




26 S7 
40 OF 




Jump to residual 


55 




L7 OF 
L2 (507) 






56 




56 (100) 
L5 OF 




Rescale 


57 




66 (507) 
S5 F 






58 


(8) 


40 OF 
L4 S5 


by 43, 74 




59 




40 IF 
19 IF 






60 




50 IF 
TO (514) 






61 


(9) 


00 IF 
LO S5 


by 45, 74 




62 


(6) 


40 S5 
L4 S5 


A>J k ; by 44, 75 

by 41, 75 




65 




40 IF 
LL IF 






64 




56 (12) 
26 (100) 




Rescale 


6 5 


(12) 


L5 (515) 

52 (15) 


from 64 




66 


(15) 


26 (14) 










50 OF 


from 65 




67 




L5 (516) 
74 OF 






68 




L4 (515) 
40 (515) 






69 




S5 F 
40 (516) 






70 


(14) 


F5 (517) 
40 (517) 


from 66 


_ 



LOCATION 

n 

72 
73 
7^ 
75 
76 
77 
78 

79 
80 
81 
82 
83 

85 
86 

87 
88 



ORDER 



NOTES 
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L8 



(11) 

(15) 



<19) 
(5) 



LO (505) 




36 (15) 




L5 (6) 




L4 (519) 




40 (6) 




F5 (8) 




42 (8) 




42 (9) 




LO (524) 




36 (18) 




00 IF 


from 71 


L5 S5 


by 44,80; from 8l 


86 iif 




00 S4 


by 40,78 


F5 GO 




40 (4) 




LO (503) 




36 (19) 




F5 (11) 




42 (11) 




LO (525) 




32 (11) 




26 (17) 




00 F 




85 11F 


from 79,88 


00 S4 


by 39, 87 


L4 S3 




40 S3 


by 41,86 


L5 (5) 




L4 (519) 




40 (5) 




F5 (19) 


.- 


40 (19) 




LO (504) 


■ 


36 (19) 




L5 (516) 





Record 128 components of An , 



j LOCATION 




ORDER 




NOTES PAGE 7 L tf 




89 




40 OF 
L5 (515) 








90 




52 (20) 
22 (24) 








91 


(20) 


00 IF 
50 (20) 


from 90 






92 




26 (Rl) 
40 OF 








95 




L5 (525) 
56 (22) 




is rt =0? 
k 




94 




L5 (525) 
LO (515) 




' 




95 




56 (25) 

L5 OF 




1 l«( r, 1 | - end constant < 0? 

1 ' k ' ' 




96 




LO..(523)' 
56 SB 




*"~C - 91 when done 1 




97 


(25) 


L5 (525) 
4o (522) 


from 95 


\ 

< 




98 


(22) 


L5 OF 

4o (525) 


from 95 


i 




99 




22 (24) 
00 F 




ll^? k+1 ll = 1 l*r k | | in ace, 




100 


(100) 


L5 (526) 
40 (101) 


from 56,64 


1 




101 


(105) 


4l OF 
50 (102) 


f rom 106 






102 


(101 ) 


7J S5 










j^O S3 


by 100, 104 






105 




L5 (101) 
L4 (519) 








104 




40 (101) 

F5 OF 




■ 




105 




kO OF 
LO (505) 




. . 




106 




56 (104) 
22 (103) 







LOCATION 




ORDER 




NOTES PAGE 8 


l8 ' 


107 


(104) 


50 (102) 
7J (506) 


from 106 






108 




ko (506) 
26 (25) 








109 


' (102 ) 


kO F 
00 F 








110 


(500) 


85 11F 
00 s4 








111 


(501) 


86 11F 
00 Sk 








112 


(502) 


85 11F 












00 F 


by 11 


85 11F 00 mk 




113 


(503) 


86 11F 












00 F 


by 10 

: 


06 11F 00 NS4 




114 


(504) 


05 11F 












00 F 


by 12 


05 11F 00 NS4 




115 


(505) 


00 F 
00 s6 




00 F 00 NF 




116 


(506) 


00 F 












00 F 


by 0/l08 


Scaling factor 




117 


(507) 


00 F 












00 F 


by 4 


a/2 + b/2 = l/tf 




118 


(508) 


00 F 
00 F 




zero 




119 


(509) 


00 F 
00 F 


by 71 X 






120 


(510) 


00 F 












00 F 


by 8 


\ 2 /2 




121 


(5H) 


80 F 
00 F 




-1 




122 


(512) 


00 F 
00 70F 








123 




00 F 












00 F 


by 26 

1 


End constant 





LOCATION 




ORDER 
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12k 


(51*) 


00 F 










00 F 


by 28,31,36 


~\ 


125 


(515) 


00 F 










00 F 


by 68 


1 1 , -> . ,2 

iN r w.iil M 


126 


(516) 


00 F 










00 F 


by 69 


I^W^L 


127 


(517) 


00 F 










00 F 


by 36,70 


Equation counter 


128 


(518) 


L4 S3 
kO S3 






129 


(519) 


00 IF 
00 IF 






130 


(520) 


00 S5 
00 S5 






131 


(521) 


S2 (7) 
40 128S5 






132 


(522) 


00 F 










00 F 


by 99 


H^k-iH 


133 


(523) 


00 F 










00 F 


by 28, 98 


ll** k ll 


134 


(524) 


NO F 
lh 128S5 






135 


(525) 


80 IF 
L5 128S5 






136 


(526) 


7J S3 
ho S3 

OOK(Rl) 







LOCATION 



ORDER 



3 
4 

5 

6 

7 



10 
11 
12 

13 
14 

15 
16 



(A) 



00 K 
50 33L 
50 L 
26 (N12) 
L5 30L 
42 4L 
42 % 
50 36OF 
50 3L 
26 (N12) 
L5 F 

86 iif 
00 2560F 

F5-5L 
40 5L 
F5 4L 

40 4L 
LO 31L 
32 4L 

41 29L 
L5 F 
l£ 29L 
40 29L 
F5 9L 
40 9L 
"L0"32"lT 

32 13L 
22 9L 
L5 5L 
40 15.L 
Fl 29L 
00 F 
00 F 
F5 5L 
40 5L 



NOTES 
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Read in a 5 "o, s r 



Read in one row of matrix and transfer 
it to drum 



Check sum 



Computation 



Form x. + S. , = S. 
1 ' 1-I 1 1 



r s_ ± = 



Store -S - 2~ 59 = CKS 
n 



on drum just after last elt. of row 



LOCATION 




ORDER 
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17 




F5 (C) 

ko (c) 






18 




LO (N) 

36 20L 




) 


19 




22 IL 
OF F 






20 




50 36IS6 
50 20L 




Read in s y and initial approximation 


21 




26 (N12) 
26 22 L 






22 




50 33 

50 22 L 






25 




26 (N12) 
L5 33L 






2k 




kO F 










L5 3^L 






25 




kO IF 
L5 35L 






26 




40 (102) 

26 20F 




Jump to iteration routine 


27 


(c) 


00 F 

00 F 






28 


(N) 


00 F 
00 s6 






29 


k 


00 F 
00 F 






30 




00 F 

00 360F 






31 




K6 (N12) 
L5 36OS6 






32 




kl 29L 
L5 36OS6 






33 




00 F 
00 F 






3^ 




00 F 
00 F 







LOCATION 
*35 



ORDER 



3 

4 

5 
6 

7 
8 

9 

10 
11 
12 

13 
14 

15 



(B) 



00 


F 


00 


F 


00 


K 


4o 46L 


K5 


F 


42 


37L 


L5 


F 


40 


(I) 


50 


(N+l) 


75 


(I) 


S5 


F 


L4 4ol 


40 


7L 


L4 


(N+l) 


40 4lL 


L5 


30(A) 


42 


8L 


00 


F 


00 


F 


22 


8L 


40 


F 


F5 


8L 


40 8L 


F5 


7L 


40 


7L 


to 4lL 


32 


12L 


26 


7L 


41 29(A) 


L5 


30(A) 


42 14L 


L7 29(A) 


L4 


F 


40 29(A) 


F5 


14L 
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L 8 



Store scale factor; 



Plant link 



Store i 



i(N+l) 2" 59 + 2560 



End constant for drum transfer loop 



Drum read order inserted by 4 



Transfer ith row of matrix, including CKS 
to W.M. 36O, 361, ... 360+N 



Form CKS 



LOCATION 



16 

17 
18 

19 
20 

21 
22 





23 




2k 




25 




26 




27 




28 




29 




30 




31 




32 




33 



ORDER 



kO l^L 
LO 42L 
36 18L 
26 l^L 
L5 8L 
LO h3L 
k2 20L 
F5 26(A) 
22 20L 
Ik Y 
kO F 

■L3-F 
36 23L 
FF F 
hi' (MSP) 

>9,'(LSPy 
L5 30(A) 
te 26L 

"15 ^8l 

Ij-6 27L 
L5 (ISP) 
50 F 

7^ F 
Ik (MSP) 
hO (M5P) 
S5 F 
*K) (LSP) 
L5 27L 
L4 A5L 
46 27L 
F5 26L 
fe 26L 
LO kjL 
32 33L 
26 26L 



—1- 



JTllUJi j- J 



Old CKB = Nev CKS? 
> yes; < No I 



Compute residual, 



double precision 



1 

LOCATION 




ORDER 
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3* 




L5 ^8L 
I* (I) 
h2 35L 




Vi = s < s J 


35 




50 ^+6l 
75 F 






36 




66 35(A) 

S5 F 






37 




LO (MSP) 

22 F 






38 


(I) 


00 F 
00 F 






39 


(N+l) 


00 F 
00 1S6 






ko 




85 11F 

00 2560F 






hi 




00 F 
00 F 






k2 




L7 29(A) 
IA 36OS6 






hi 




00 IF 
00 IF 






kk 


(M3P) 


00 F 
00 F 






*5 


(ISP) 


00 F 
00 F 






k6 




00 F 
00 F 




s, stored here 
k 


hi 




L5 (ISP) 
50 36OS6 




End constant 


k8 




00 S3 

00 361S6 







LOCATION 




ORDER 
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00 K 











L5 (506) 
26 1L 






1 




50 12F 
50 1L 




Print out final scale factor 


2 




26 (Pl6) 

L5 5L 






3 




I* (n) 

kO 11L 






k 




92 135F 
92 515F 






5 




22 5L 
L5 S3 






6 




50 S9 
50 6L 






7 




26 (Pl6) 
92 131F 






8 




92 515F 
F5 5L 






9 




kO 5L 
LO 11L 






10 




32 19(A) 
22 5L 






11 




00 F 
00 F 

(Pl6) OOK 

(N12) OOK 

W(k) 26lN 

















L 8 



